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NATIONAL ADVISORY COMMITTEE FOR AERONAUTIUS ` 
TECHNICAL NOTE 3288 


ON THE ANALYSIS OF LINEAR AND NONLINEAR DYNAMICAL SYSTEMS 
FROM TRANSIENT-RESPONSE DATA 


By Marvin Shinbrot 
SUMMARY 


A general theory of the so-called "eguations-of-motion” methods for 
the analysis of linear dynamical systems is developed first. It is then 
shown that when viewed from this general point of vantage, all of these 
linear methods can be extended in a straightforward manner to apply to 
the analysis of nonlinear systems. In addition, through use of this 
theory, & new method is derived. It is essentially & variation of the 
well-known "Fourier transform" method for the analysis of linear systems 
but possesses certain advantages over previous methods. Application and 
effectiveness of this method are demonstrated by three examples, two of 
which are nonlinear - one highly so - and the third being of the fourth 
order. 


INTRODUCTION 


It has often been suggested (e.g., in ref. 1) that nonlinearities 
which are ignored in the classical theory of the equations of motion of 
an aircraft may be responsible for certain unusual phenomena which have 
been observed in flights of modern high-speed airplanes and missiles. 
Consequently, it seems desirable to develop methods for the analysis of 
such nonlinear systems - methods which allow the calculation from measured 
transient-response data of the nonlinear stability characteristics as well 
as the classical linear stability derivatives of the aircraft. Several 
such methods are described in reference 2, the principal one consisting 
of a generalization of the so-called “derivative method" which was orig- 
inally devised for use with linear systems (cf. ref. 3). However, the 
methods described in reference 2 leave something to be desired from both 
points of view of accuracy of the results and lengthiness of the calcu- 
lations. In addition, application of these methods requires, in all but 
the simplest cases, the previous evaluation by some means of those sta- 
bility characteristics which are linear. In view of these shortcomings, 
an attempt has been made in the present study to find simpler, more accu- 
rate, and more general procedures. The problem is attacked by first exam- 
ining several well-known methods for the analysis of simple linear systems 
and then modifying them 88 necessary to allow their application to more 
general systems. 
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Many methods for the analysis of linear systems have been proposed 
in the past (see, e.g., refs. 3 and 4). In reference 5, these methods 
have been classified under two main heads:  l"equations-of-motion" methods 
and "response-curve-fitting" methods, the former title including the 
derivative method and what have been called the Laplace transform and the 
Fourier transform methods (ref. 3), and the latter consisting of such 
methods as Prony's (refs. 3 and 6) and the techniques of reference ἢ, 
Since the response-curve-fitting methods involve the explicit solution 
of the equations of motion in terms of the physical parameters of the 
system at hand, they do not seem suitable for use with nonlinear systems. 
Hence, we sh&ll.be concerned solely with the equations-of-motion methods. 


Each of these methods has been considered in the literature as an 
independent. entity; apparently, no attempt has ever been made to subsume 
all of them under a single general theory. For the purposes of the pre- 
sent study, such a theory would be desirable since it seems reasonable to 
expect first that when viewed from a more general point of view, a gen- 
eralization of the methods to nonlinear systems might appear; and second 
that once such a theory is known, it might be possible to develop new 
methods, superior in certain respects to the old ones from which the 


theory sprang. 


In accordance with this pian, the paper begins with a short presen- 
tation of the three best known of the equations-of-motion methods. These 
methods are examined from a new point of view which is then shown to lead 
to the general theory for linear systems. The further extension to non- 
linear systems is considered next. Based on the general theory, develop- 
ment of a new method for data analysis follows. Finally, some examples 
of the application of this recommended method are given. 


SYMBOLS 
ον angle of attack, radians Е 
š да | 
at { 
> p-k 
i y y 
Š Mẹ 15М4 
O  —r цал «μα 
c Lo 
1 ~ mV 
Š elevator deflection, radians 


Ту pitching moment of inertia, slug-f+2 
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k 2 Ма, πα М 
L linear lift force, lb 
aL 
5 δα 
OL 
- X 
L(a) nonlinear lift force, lb 
m mass of aircraft, slugs 
M linear pitching moment, ft-1b 
ОМ 
ын δα 
. ОМ 
- δά 
Ma SM 
95 
ОМ 
Ма dd 
M( 0) nonlinear pitching moment, ft-lb 
q pitching velocity, radians /sec 
t time, sec 
V velocity of aircraft, ft/sec 


Other symbols will be defined as they are introduced, 


ANALYSIS 


Some Equations-of-Motion Methods for Linear Systems 


in the presentation of the general theory of equations-of-motion 
methods, it is desirable to have some examples of these methods set down 
for examination. The three best known of these methods are briefly 
described below; for a more detailed discussion of them, see reference 3. 
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Ав a concrete example, let us consider ап airplane operating under 
conditions where the stability characteristics are effectively linear, 
so that, as in reference 22% the di d its longitudinal ΕΕ can 
be written | | 


- aly - ашу + quí = Dig 
- OM, - ἀΜὰ - gMq + QIy = ὃΜε 


It will be convenient to eliminate q and so to write these equations as 
alt) + ba(t) + kx(t) = Colt) + С,5(%) (2) 


where dots denote differentiation with respect to time +, It is assumed 
that time histories of a(t) and 5(t) are available from flight records 
and that one wishes to calculate the constants b, k, Co, and σι. For 
simplicity, it will further be assumed that the μι... b is "positive 
and that (t) is a pulse, во that (+) is zero for t sufficiently 
large. These restrictions will be removed farther on. 


The derivative method.- In order to apply the derivative method, it 
is necessary first to differentiate the given records to obtain the 
derivatives a(t), a(t), and 5(t). Then, for fixed t, equation (1) may 
be considered as an equation for the constants b, k, Со» and Cy. By 
letting t vary, a set of such equations can be obtained, Which can in 
turn be solved by least squares (see reference 6, page 210, where, however, 
the integral rather than the sum of the squares of the errors has been 
minimized) for the desired constants. By this procedure, the following 
equations are obtained: 


со 
> / авав + / “аав - со [ авав - o, / ἆ bare -/ ὁ ἃδε 
Ο 
со со со o. co 
of aadt + κ’ оға - co / adt - οι / aodt = -/ адаб 
O ο ο О О 
со со со со : oo 
- of Sadt - x f Sadt + co | 52а + o f 5946 -/ Sadt 
ο ο O о ο 
° o. o. 90, =, 
of Š dat -x f боле + Co f бба + e f 2а -/ Š adt 
ο ο 


О O ο 


(9) 


Equations (2) can be solved for the desired parameters b, k, Co, and C,. 
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The Laplace transform metkod.- Letting A(p) and Alp) denote the 


Laplace transforms of a(t) and 5(t), respectively, во that 


A(p) = «ρε o( t)dt 
л ” 


А(р) - е-Р© 5(+)а+ 


ο 
it follows that if о(%) and 5(t) are related by equation (1), then 


(р + bp + k) A(p) = (Cyp + Co) Alp) (4) 


(ref. T). In writing down equation (4), it has been assumed for simplic- 
ity that о(0) = à(0) = 8(0) = О; this restriction is inessential and will 
be removed later on. For any value of p, equation (4) is an equation 

in b, k, Co, and 61. After finding the Laplace transforms of a(t) and 
S(t) for several such values of р, say for р =P, Pas + + ., Py, the 
corresponding equations (4) can be set up and solved by least squares to 
obtain ee 


b y py2A2 (p) +k ). p,A7(p4) - Со 1 pyA(p3) А(р;) - 
Са ) ру?А(р;) A(p1) = - ) pí А (py) 
b у pi (p + k ) βίο) -Co ) Ам Абы) - 
с. ) paA(pi) δι) = - У arar to) 
-b Уу νι) дн) - к ) Αρ) Ат) + Co ) (pa) + 
с, y р,42( ғ) = у p42A(04) д(р;) 
- b уығы) Alpi) - Е > p4A(p4) A(p1) +Co Y p4A (p1) + 


с, ) ρα δ (pi) = ) ρα A(pi) Alpi) 


(5) 


where all sums are over the range 1, . . ., N of the index i. These 
equations then can be solved for the desired constants. 


The Fourier transform method.- The Fourier transform method proceeds 
in much the same way as did the Laplace transform method. 
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Defining | | dco 5 
A(iw) - Г sum alt) dt 
ven (6) 
Aliw) = | ο-λωῦ s(t) at 
J, 
(12 = - 1), it follows that-if α(ο) = à(0) = (0), then 
[(10)% + iub + k] Α(1ω) = [0ο + 1061 Aliw) (7) 


This equation can be written as two real equations by setting 


A(io) = σ(ω) - iS(w) 
(8) 


Aliw) = Γ(ω) - i E(w) 


Putting equations (6) and (8) together, we see this means that 


σ(ω) - [Г оа) cos wtat S(w) = fal) sin wtdt— 


о 


со 
A (t) sin wtdt ` 


O 


Γ(ω) - sc) cos wtdt Σ(ω) 


Breaking equation (7) into its real and imaginary parts gives 


w S(w)b + C(w)k - Mw)Co - e E(w) С, = Cw) 


w σ(ω)ο - S(wW)k + X (w)Co - w T(w)c, = - 5 8(0) 


After C(w), S(w), Г(ш), and Σ(ω) have been evaluated for several values 
of w, substitution into these equations yields a number of equations 

for b, k, Co, and C, which can be solved by least squares. These equa- 
tions, corresponding to (2) and (5), are exceedingly complicated and will 
not be reproduced here. жин | | Е 1 
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The Common Feature of These Methods and the First Generalization 


Each of these methods is usually derived as in the preceding section, 
through use of a certain specialized concept. Thus, the derivative method 
18 based on the fact that for fixed +t, equation (1) is linear in b, k, 
Co, and C,, so that a set of simultaneous linear equations for these 
parameters can be obtained by varying t. The Laplace and Fourier trans- 
form methods stem from the theory of these operators. Thus, if a(t) 
and 8(t) are related by the linear differential equation (1), their 
Laplace transforms are related by equation (4). However, equation (19) 
is linear in the parameters, and so by writing this equation for several 
values of p, one can solve the resulting equations by least squares for 
the desired constants, These derivations tend to obscure the common 
idea which can be shown to lie behind all the methods. This difficulty 
can be overcome if these particular derivations are forgotten and if 
attention is fixed entirely on the formal processes whereby the final 
least squares equations are obtained. With this in mind, let us recon- 
sider the three methods. 


The derivative method.- Equations (2) for the derivative method are 
formally obtained by 


(1) Multiplying equation (1) by the four functions 
a(t), a(t), S(t), and (t) one at a time and 
(2) Integrating the results from zero to infinity. 


We should forget for the moment the interpretation of this procedure as 
the solution of equation (1) by least squares, and simply keep the process 
of multiplication and integration in mind. 


The Laplace transform method.- A similar process can be described 
for the Laplace transform method. Choosing N(>4) positive numbers pj, 


(1) Equation (1) can be multiplied by the functions 


е-ріб, і = 1, . . ., N, and 
(2) The results can be integrated from zero to infinity. 


If the resulting equations are solved by least squares, precisely equa- 
tions (5) for the determination of the parameters by the Laplace trans- 
form method are obtained, provided that in step (2) any integrals which 
arise involving derivatives of a(t) and 5(t) are integrated by parts 
to eliminate these derivatives. 


It should be noted that although there appears to be one more step 
here than there was in applying the derivative method - notably an addi- 
tional least squares following step (2) - this addition is more apparent 
than real, since it is necessary to apply a least squares process here 
merely because М (which is generally greater than four) equations are 
obtained for the four parameters, while in the derivative method exactly 
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four such equations were obtained. Thus, the least squares step could 
be eliminated by choosing N = Y (of course, such a choice is not really 
practical), or, alternatively, such a step could be added to the deriv- 
ative method by solving the equations resulting from step (2) in that 
method by least squares instead of in the ordinary way. Of course, such 
a step would only be made to make the two methods so far described 
formally more similar; in practice, it would not be performed. 


The Fourier transform method.- Finally, 


(1) If equation (1) is multiplied by cos wt and sin wt 
for several values of w, and 

(2) If the results are integrated from zero to infinity 
(as in the Laplace transform method, integrating by parts 
to — explicit dependence on the derivatives of q 
and 6), 


one obtains a set of equations identical with those obtained from the 
Fourier transform method. - 


The general method for linear systems.- The general development of 
equations-of-motion methods is now manifest. One takes the equations of 
motion for the physical system under consideration - for definiteness, 
say equation (1) - and 


(1) Multiplies them by N arbitrary (but sufficiently 
smooth) functions yy(t). 

(2) The resulting equations are then integrated between 
two definite limits, say, zero and Т, 


In the three methods described above, T = о, but this is not essential. 
In order to avoid some complications initially, we shall continue to 
integrate over this infinite interval; this restriction will subsequently 
be removed, however, and T will be allowed to have finite values. In 
the case of equation (1), the process just described leads to N equa- 
tions of the form 


f yy (tJa(t)jat + xf yy(t)e(t)at - cof yy(t)8(t)at— 


О О O 


Ci [ yy (5)8(t)at-- 47” yy (t)&(t)at— V21,...,N (9) 


O O 


It is possible that the functions yy(t) depend on a(t) or S(t) as, for 
example, in the derivative method; in such cases, equations (9) can be 
considered as N equations which are to be solved by least squares for 
the desired parameters. Of course, this process requires the calculation 
of the derivatives a(t), &(t), and S(t). On the other hand, if the . 
functions yy Ct) are explicitly independent of a and 5, as is the case 
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in the Laplace and Fourier transform methods, the following formulas, 
obtained by integrating by parts, are used: 


f yy(tja(t)dat 


О 


Ј эма Da = - уу(о)&(о) + ὑψίο)αίο) + [7 Ууссан 


O O 


J жайға = - умож(о) - / ыда 


= vylo}a(o) - f^ ὑψίύλαίολάι 


Substitution into equation (9) gives 


: ъ| υγ(ολα(ο) ε΄ jy Cos Coo | +x yq) - 


O O 


co Υνίθδ(θαν + o, уу бово) + ὑνίε)οίο)ας | = vylodilo) - 


y y(Oja (0) n yy(tja(tjat, v = 1, 2, . .. +s N (10) 
О 


Equations (10) аге N equations in b, k, Co, and Ci. If N24, they 
m&y be solved by least squares for these parameters. 


The choice of the functions y(t) defines which equations-of-motion 
method is being used. For this reason, these functions will be referred 
to as the "method functions.” 


Generalization to Nonlinear Systems 


In this section, all considerations will refer to the equations of 
longitudinal motion of an aircraft. It will be seen that the method 
actually is applicable to a far wider class of equations - in particular 
to the equations of lateral motion of an aircraft, including, if it is 
so desired, cross-coupling terms. The special analysis presented here 
can be extended to other problems, the only real restriction being that 
for practical reasons too many parameters cannot be handled at once. 


The following equations, which involve assumptions of constant air- 
Speed, smallness of certain quantities, etc., &re often used to describe 
the motions of an aircraft which has linear stability characteristics: 


- oly - amV + qmV = δίς 
- OMe, - GM, - qMq + GIy = δΜς 
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If, on the other hand, it is assumed that the 117% and moment functions 
are nonlinear in a, these equations become 


- По) - дау + quí = διε 
(11) 
- Μία) - амд - «М0 + Шу = 5М 


We shall assume that over the range of interest, approximations of the 
following form are valid: 


(12) 


Μ(ω) = Μια + Маё +. . . + Май 


il 


where the coefficients 14, М; are constant, Only the first three terms 
of this series will be retained in the present analysis since this three- 
term approximation usually balances very well the opposing requirements 
of simplicity and adequate representation of the aerodynamic parameters. 
If more terms are found to be necessary in a particular problem they can, 
Of course, be added. It should be noted that this three-term арргохїша.- 
tion 18 still fairly general, even кетише the possibility of ΕΠ 
in the nonlinearities. 


By use of the approximations (12) with 1 = n = 3, equations (11) 
can be written 

Li La 4 Ls з i1 Le 

"av 7a “шу? ce tuy 

(13) 
M M А . 

- 22 9 τα са „Ма ав. ta - hat Gnas 

y y y y y y 


The generalization of the methods of the preceding section to such non- 
linear systems proceeds in the obvious way. First, multiply each of 
equations (13) by N method functions y,(t) which have been selected 

as suitable. This operation is followed by integration of the resulting 
ОМ equations. If records of the derivatives qanda@ are not available, 
one then eliminates the terms involving these derivatives by integration 
by parts. The result is М equations in Τη» L5, Lg, and Ig and N equa- 
tions in M,, Ma, Mg, Уі, and Mg. If Ν2 5, these two sets of equations 
can be solved by Least squares for the parameters. 


Of course, the pitching velocity q can be eliminated from equa- 
tions (11} to yield the single equation ян 
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Mi M Ea) pa _ Ma) L(a) _ M5 + Ma le 
цан | тоо τν | ш Iy o> mV T Ty mV 
Lt(a) = em = (а) 


With the approximations (12), this leads to the following generalization 
of equation (1): 


& + (bo + bja + Бао )6 + (ko + k,0 + ksa%)a = CoS + C,5 (14) 
where 
bo = MG, qd ae a ko = (E + Mg La 
_ 212 | _ Ма, Μα le 
bi = aF dg Vi = 025) 
y iy 
3Ls _ » Mg Із Ls 
ра = шү SET το» Ty αὐ 


By applying the method described to equation (14), the constants bo» Ві, 
Ыр, Ко, ky, Ка can be calculated, 


Solution for the parameters in equations of the general form of 
equation (14) is of interest to workers in many fields. In addition, 
the evaluation, from these parameters, of the stability constants occur- 
ring on the right sides of equations (15) is of considerable interest to 
aeronautical engineers, and so this problem will be considered in further 
detail. One cannot, in general, isolate the constants 14 and Mj to 
obtain, from values of the bi and the ky, the nonlinear functions L(a) 
and Μία). For this reason, it is ordinarily best to apply the method 
directly to equations (13) rather than to equation (14), provided that 
records of both q(t) and a(t) are available. 


On the other hand, there is one case of interest when equation (15) 
can be used directly and measurements of q(t) are not needed. It some- 
times happens that while the pitching moment Μία) is nonlinear » the lift 
L(a) can still be successfully approximated by a linear function. In 
this case, we have the approximations 


L(a) = Lye 


(16) 
М(о) = Μια + Ma? + Mago? 
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and so equation (14) becomes 


б + b& + (Ko + Έα + κραξ)α = CoS + Суб (17) 
where 
€ Μα La 
add бге 
u да 
"ЭРЭР (18) 
ly 
Ms 
Κα = = =. 
Ly 
“Ма la | 
Thus in this case, M(a) can, except for the term == ny occurring in the 


expression for Ко, be obtained from an analysis oF equation (17). At 
high speeds, however, this term is small and its effect on the curve of— 
Μία) versus a can be neglected. Thus, the expression for kj in 
formulas (18) can be replaced by the expression Š 


M 
πμ. (19) 
y 


and in this case, the nonlinear moment M(a) is completely determined by 
the knowledge of Ко, Κι, Қа, and, of course, Iy. 


Choice of the Method.Functions 


Up to this point in the. general discussion, the method functions 
yy (t) have been to a very great extent arbitrary, having to satisfy only 
certain weak smoothness conditions. In this section, the possiblity of 
developing new and perhaps. improved methods for data analysis by means 
of a particular choice of the method functions will be explored. 


Previous experience, consisting in part of unpublished analyses per- 
formed at Ames Aeronautical Laboratory, have indicated that the Fourier 
transform method generally results in greater accuracy than either the 
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derivative or the Laplace transform methods.1 For this reason, the 
method to be discussed will be a variation of the Fourier transform 
method. | 


Let us consider first some evident shortcomings in the Fourier 
transform method. These defects will offer definite goals to be held 
in mind in the development of a new method. First, there is a weakness 
in the Fourier transform method in that all integrations proceed over 
the interval from zero to infinity (cf. eq. (10)). This causes diffi- 
culties in any example in which a(t) and 5(t) do not approach zero so 
rapidly that their integrals exist. Thus, referring, for example, to 
equation (1), if Ὁ is not positive or if 5(t) does not approach zero 
quickly enough, the method cannot be applied straightforwardly. Further- 
more, even if Ὁ >0 and 8(t) ->0 іп such a way that fa(t)at exists, 
the experimental record often is not long enough for this integral to be 
accurately calculable when the system is so lightly damped (b smell) 
that sizable oscillations persist even to the end of the run. One device 
which is sometimes used to overcome this difficulty is equivalent to a 
change in the method functions. Instead of the functions sin wt and 
сов Ut, the functions e-B$ sin wt and e"Ët cos wt, with some fixed 
constant B, are used. However, this leads to the same objection that 
was voiced in footnote 1 for the derivative and Laplace transform 
methods, notably that the method functions approach zero. Other tricks 
for dealing with such deficiencies in the Fourier transform method can 
be evolved, but, rather than develop new devices for each special case, 
it appears wiser to construct a generally applicable method in which 
these difficulties never arise - that is, one in which the integration 
proceeds only over a finite interval. 


The second defect which we shall consider becomes clear from an 
inspection of equation (10), with the functions yy(t) of the form 
sin wt and сов wt for certain values of ш. Referring to equation (10), 
it is easily seen that one point, namely the point t = 0, is weighted 
very heavily because of the occurrence of the quantities a(0), a(o), 
and 3(0). It should be noted that not only are the values of a(t) and 
S(t) at one point relied on to this great extent, but that even the 
relatively inaccurate value of the derivative of a(t) at that point is 
weighted. Thus, advantages in accuracy might be expected to accrue if 
these terms were eliminated. 


Since the method for overcoming this second deficiency in the Fourier 
transform method will also be used in treating the first, the problem of 
eliminating dependence on the initial values will be discussed now. To 
this end, consider equation {10). We begin with the Fourier transform 


ltt would seem that no rational explanation for this conclusion has 
heretofore been offered. However, the theory described herein appears to 
afford such an explanation. A long and rather tedious analysis based on 
this theory has indicated that the failure of both the derivative and the 
Laplace transform methods is due in large part to the fact that the asso- 
ciated method functions approach zero very rapidly as time progresses. 
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method, that is, with method functions yy(t) of the forms sin wt 

and cos wt. Noting that most-of the terms which depend on the initial 
values at t= O are multiplied by уу(0), it is seen that a choice of 
method functions such that уу(0) is zero for all у = 1, 2, . . ., N 
represents a step in the right direction. Such a choice is easy to 
make, since it is only necessary to eliminate those method functions 
which have the form cos wt; then, we may write 


This does not entirely eliminate the dependence on the initial conditions, 
however, as the term Уұ(О)о(0) remains in equation (10). If this term 
can also be removed, the second weakness in the Fourier transform method 
Will have been entirely corrected. This will clearly be the case if 
ὑγίο) as well as уу(0) is zero forall У-1,..., N. A possible 
choice of the method functions for which this is so, a choice which 

still retains the advantages of the favored Fourier transform method, is 
the following: | 


1 - сов 20у 
yy(t) = sin?wyt = — qom ies Νὶ (21а) 


With this choice of the method functions, equation (10) becomes 


- o f a(t)yy(t)at + к [7 a(t)yy(t)at - co f S(t)y (t)dt + 
O O ο (22) 


e f. S(t)yy(t)at = - | at), (tat 


о о 


The method functions (21а) would be used for systems satisfying dif- 
ferential equations, like (1), which involve derivatives of the second 
order. More generally, and for the same reasons, if the highest order 
derivative occurring in amy of the equations of motion of a system is 
the nth, the following method functions are suggested: 


yy(t) = sin? wyt, V = 1, +< + ə М (21b) 


Thus, in the case of the two-degrees-of-freedom system described by 
equations (13), there is no point in using formula (21a); the simpler 
method functions (20) may as well be used. 


As for the first of the weaknesses in the Fourier transform method, 
that which is due to integration over an infinite interval, it would 


NACA TN 3288 19 


appear at first glance to be easily disposed of by merely choosing some 
finite, positive number T and integrating over the interval from zero 
to T, This, however, introduces another difficulty. Suppose yy(t) is 
given by equation (21a), so that yy(0) = yy(0) = O. Returning to the 
derivation of equation (10), it may be seen that if one integrates only 
over the interval O<t<T, all the good which has been achieved by 
eliminating dependence on the point t = O is obviated by certain terms 
which arise - terms of the form yy(T)a(T), уу(Т)а(Т), and yy(T)e(T). 
Thus, the heavy dependence on the initial conditions is replaced by a 
dependence on the final conditions. Of course, the same approach as 

was used for eliminating the initial conditions can be used again - that 
is, the method functions can be chosen in such a way that 

yy(T) = y (T) = O. One possibility, naturally, is to choose the frequen- 
eles wy such that 


sin юу? = O, М-і,..., М 


Thus, we cam set 


This choice of the frequencies (corresponding to the method functions: 
given by equation (21b)) leads to an elegant method which gives satis- 
factory results in certain cases. On the other hand, the difference 


m between two successive frequencies is too large to define the 


"Frequency response" (to use loosely the terminology of the Fourier 
transform) of some examples adequately. For this reason, we should like 
to be able to choose the frequencies as follows: 


Wy = M, v = 1, e е οὐ N (23) 


Of course, this means that for frequencies having an odd subscript, yy(T) 
will be different from zero. To overcome this difficulty, the following 
choice of the method functions can be made: If the highest derivative 
occurring in the equations of motion is the nth, define the method func- 
tions by the formulas 


You C5) = sin? ш. 


эр. 
2н 
PEE e. (24) 
sin “әд © gon ο Suri T 
Угра- (9) z 
2р 


°, 2р+1 


16 : NACA TN 3288 


where the frequencies Oy are given by equation (23) and Т is the 
length of the run. What has been done in choosing these method functions 
is the following: Those method functions which are such that yy(T) is 
zero (1.е., those method functions which have even subscripts), have been 
left unaltered in the form given by equation (21b). The remaining method 
functions have had the last quarter cycle, during which they would nor- 
mally have varied from O to +1, chopped off, so that they are identically 
zero over p&rt-of the interval О < t < T. The claim is not made that- 
this is the best of &ll possible choices for the method functions. There 
18 certainly no reason for such an &ssertion, particularly in view of the 
fact that & certain amount of the data is not used by each of the odd- 
numbered method functions. However, this amount is small after all and 
no d&tum is completely discarded, since the even-numbered method functions 
use all the data, That these method functions do seem to be adequate is 
indicated by the results obtained in the examples given below. 


Before-proceeding to the examples, one further change, imposed in 
the interest of simplicity in the computations, will be made in the method 
functions. For the odd-numbered method functions, certain complications 
arise in the computations due to the fact that the point (82р/2р-1)Г at 
which the function is cut-off may not coincide with a point at which the 
data are tabulated. For this reason, values of even-numbered frequencies 
will be chosen in accordance with equation (23). The odd-numbered fre- 
quencies, on the other hand, will be changed (by as little as possible, 
to be sure) in such a way that the following condition is satisfied. Let 
the data be tabulated at the points t = to; t1, + + ., tg, where 
tin, - tk = constant. The odd-numbered frequencies are then assumed to 
be chosen such that | 


sin ὠρμρτῖα = O 


where tx is that one of the two tabulation instants closest to the time 
(2u/2u+1)T having an even subscript. The-reason for this last condition 
is merely that it 1s convenient for a numerical integration procedure 
(such as the one in the Appendix) to have an even number of intervals 
(tk, ар); thus, Simpson's rule, for example, calls for an even number 
of intervals in its application. 


For many of the experiments performed on airplanes and missiles, the 
run is 2 seconds long and the time between tabulated points is 0.05 sec- 
ond, go that the tabulation times t, are at О, 0.05, 0.10, . . ., 2.00 
seconds. For these values of Т and At, the rule given above for the 


frequencies ων ів climaxed by the following table: 
ΕΠΕῚ {51611189 [10| 11 [12| 13 | as | 15 
ων Цас 


τ || || ое (сол | ae | con 
2113 Tl 2 | 17 9 } 5 | 9 3 | 2 19 
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The frequency w, has been omitted from the table since, according to 
the rule for the determination of the method functions, y,(t) = O 

(see eq. (24), ff.) and offers no information. A 2-second interval will 
be taken as standard in this report, and all computations will be based 
on such an interval. If one has a data run T seconds in length, it 
cam of course be brought into this standard form by a preliminary sub- 
stitution of the form 


and we always assume such a transformation has been made. 


Thus, finally, the method may be summarized as follows: Select N 
frequencies in accordance with the rule given in the preceding paragraph. 
(The number N is chosen, as in the Fourier transform method, large 
enough to cover the frequency range of interest in the particular problem; 
usually, N = 16 is adequate.) Multiply the equations of motion by the 
method functions (24) and integrate the resulting equations from zero 
to T, where T denotes the length of the data run. Eliminating all 
explicit dependence of these equations on derivatives of the data by 
successive integrations by parts results in N linear simultaneous egua- 
tions for the parameters. The coefficients in these equations are all 
integrals involving the recorded data; after these have been evaluated 
by some means,“ the equations can Бе solved by least squares for the 
desired parameters. 


EXAMPLES 


Three example problems will be solved in order to demonstrate the 
effectiveness of the proposed analysis method and to illustrate asso- 
ciated computing techniques. An effort was made to select examples 
representative of problems which often occur in aircraft-response flight 
testing and which have not been handied adequately by other known analy- 
sis methods. These examples have been simplified in some respects, not 
because of fundamental limitations of the method, but in order to avoid 
obscuring the essentials of the method and of the related computing 
techniques. 


Although limited use was made of automatic digital computing machin- 
ery in the following analyses, it did not appear worthwhile to mechanize 
complete calculation procedures for these isolated illustrative examples. 
However, the method appears to be well suited to such mechanization. 


See, for example, the Appendix, where a technique well suited for 
the type of integrations needed for this method is described. 
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Example Т 


The first example concerns the longitudinal response of a hypothet- 
ical missile for which it is assumed that the lift varies linearly with 
angle of attack, although the pitching moment does not. The velocity of 
the missile is assumed to be sufficiently high for_the expression (18) 
for kg to be simplified to i 


so that equation (17) cam be used to determine the lift and moment char- 
acteristics from a transient response. A pulse response of a system 
described by equation: (17) was. obteined from a Reeves Electronic Analogue 
Computer and it was decided to consider this response as given data to be 
analyzed. The moment.of inertia of the missile was chosen to be 100 slug- 
feet=. The nonlinear moment curve Μία) versus a which was used to 
obtain the data is shown as the solid curve in figure 1. The linear sta- 
bility derivatives were chosen in such a way that the damping parameter 

is given by | 


b = 2 - (25) 


Tm order to simplify the presentation, it was decided that free 
oscillations alone would be analyzed to determine only the constants 
occurring on the left-hand side of equation (17). Thus, for the data 
which will be analyzed, 5(t) = O, and equation (17) becomes 


o bå + (Ko + ka + қо2)а = 0 (26) 


A plot of the a(t) "data" is given in figure 2, and this information ів 
listed in table I. | 2 ы 4 2 


Since equation (26) is of the second order, we shall choose, accord- 
ing to the rule given earlier, the method functions (24) with n = 2; 
Bince the run is 2 seconds long and the time interval between data points 
is 0.05 second, the frequencies Wy are chosen as in the table on 
page 16. мэ i : - 

Integrating factors Га(уу) corresponding to each function yy(t) 
апа its first two derivatives are tabulated in columns 6 through 50 of 
table I. (As discussed in the Appendix, the UI, (yy) are numbers chosen 
such that for any integrable function x(t), the sum 


), x(n) ayy) 


n 
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is an approximation to the integral of x(t)yy(t).) Accordingly, the 
sums displayed beneath table I are the integrals needed for the reduc- 
tion of the data, divided by the factor At = 0.05. Since this factor 
occurs homogeneously throughout the 50801 from which the parameters 
are to be found (i.e., the generalization of eqs. (22) to the nonlinear 
eq. at hand), it can be divided out of these equations, and the sums can 
be used directly without first multiplying them by At. 


In accordance with the method as it has been described, the equa- 
tions which are to be solved by least squares have the form 


2 2 2 
78 EM a(t )yy(t)dt + ko | a(t)yy(t)at + 4 o? (t)yv(t)Jat + 


2 E" | 
=“ | (уо) | =- Жж αἱ τ)ϑγίο)άτ 


The sums below table I are needed for the evaluation of the coefficients 
of b, Ko, k,, and K, in the above equation. These sums are again 

listed in table II, and the coefficients in the last equation are set 

down as columns L, 5, 6, T, and 9 of table II. The sums displayed beneath 
table TT are needed for the final least squares step of the solution. 
Using these sums, it can be seen that the following equations are to be 
solved for the parameters: 


24,1669 b - 1.05270 ko - 0.336760 k, - 0.0140294 Е = - 7.01790 
- 1.05270 b + 0.885340 ko + 0.181127 k, + 0.00906192 ko = 45.5226 
- 0.336160 b + 0.144127 ko + 0.0374138 к, + 0.00189656 k, = 7.00211 
- 0.0L40294 b + 0.00906192 ko + 0.00189656 k, + 0.000108975 kə = 0.459829 


(27) 


li may seem odd to some that four significant figures have been used in 
table II for the values of the integrals and six significant figures for 
the coefficients in equations (27), while the test data are not given to 
more than three significant figures. The reason for carrying more sig- 
nificant figures in the computations than there are in the data is to 
avoid eventual loss of accuracy due to round-off errors and other errors 
of a similar type. This procedure of carrying a few more (fictitious) 
figures than the data supply is usually necessary in order to retain even 
the basic information which is in the data. 
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Solving equations (27) gives 


b = 1.95 
κ т- 50.4 

> (28) 
k, = - 30.5 

κα = 806 


The only parameter whose numerical value is given and which can be 
immediately checked is the damping parameter. Comparing the values 
of b given in equations (25) and (28), we see that 10 has been found 
with an error of 2.5 percent. The constants Ко, К,, Ка cannot be 
checked directly; however, the calculated pitching-moment curve 


fi 


Μία) = - Τγίκοα + kja + kpa?) 


- 5040 a + 3050 a? - 80600 a? 


can be plotted and compared with the true curve from which we started. 
This has been done in figure 1 from which it can be seen that the error 
at the least accurate point is less than 3 percent. 


It should be noted that the values of In(yy) given in table I can 
be used to solve any problem of the type considered here which depends 
on a second-order. differential equation or on a system of such equations. 
If the data run is 2 seconds long, it is only necessary to insert the 
data in.table I in place of the data used in this example_and proceed 
as we have just done. Ав mentioned earlier, if the data run is more or 
less than 2 seconds long, it is only necessary to make a preliminary 
transformation of the time scale so that in the new time scale the data 
run is 2 seconds long, a process illustrated in example II. 


Example 11 


The first example served to illustrate the application of the method 
to an'equation of the form (14), corresponding in the missile pitch- 
response problem to the case where only aft) and 8(t) are measured. А 
problem involving equations, like (13), of the first order, corresponding 
to the case where q(t) is available in addition to a(t) and 8(t), will 
be illustrated now. 
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Ав in example I, the lift force will be assumed linear and the 
pitching moment nonlinear. The following parametric values were 
assumed: 


m = 2 Ig, = 1000 
V=750 Ἢ Mg, = - 200 (29) 
L, = 100 Mq = - 500 


The nonlinear Μία) is plotted as the solid curve in figure 3. It should 
be noted that in contrast to the first example, the pitching moment is 
unstable at a= O and highly nonlinear. 


The "test data” were manufactured by determining a pulse response of 
this missile on the REAC. Again, the control characteristics of the sys- 
tem will not be considered, so that only free oscillations are shown in 
figure Mh. 


Merely to have a standard length of run, a 2-secqnd interval was 
always selected for the calculation of the integrating factors T (уу) 
(see the Appendix for the definition of these quantities). To illustrate 
the computation procedure for data runs of different lengths, a 1-second 
run will be considered in the present example. 


In order to use the integrating factors displayed in table IIT, it 
is necessary to make a preliminary transformation of the form 


Li le „2 _ Les da - l6 
“пу * 7 mV - у 2 dr + 2 * av 5 
(30) 
«δια. Er Boas. ϱ Wa да _ 34 Ms κ 


Recalling that for the problem under discussion, Ly = lg, Le = Ls = O, 
5(t) = O, it сап be seen from equations (30) that the equations to be 
solved by least squares for the parameters have the forms 
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2 | 2 
2 = «(т)уу(т)ат + 2 га al τ) sus dT +f а(т)уу(т)4т = 0 


O O 


2 2 2 
- ын Ñ α(τ)γγ(τ)ᾶτ - τα / a? (τ)γν(τ)ᾶτ - aj a9(T)yy(7)àT + 


мә Г? dyy(T) S 2 dyy (Tr) 
2E/ «ner. τὰ а(©зу(®ат- -/ a(t) ERC ar = ο 


| The "data" of figure 4 are presented as functions of т in 

table III. The sums, which when multiplied by AT = 0.05 approximate 
the integrals in the last two equations, are given below table III. 
These sums have been listed in the appropriate places in table IV. With 
circled numbers referring to colums in table IV, it can be seen that 
the above two equations are equivalent to the following: 


18 
© 
© 
© 


- 0 -т © -т 9+2 0-10-20 -0 
Hence, : = 
_ „22 Q x © +» @ x @ 
х @ 
or 5 
Ig, = 1033 | | üI 


using the given values of m and V (eqs. (29)), while the equations for 
the moment parameters are 


=) o CP 12:22 ) 9х9-02)) Oxo: бте 
DKE Әз) 0%+%) Ox@- Св) әне еке» 

-2) @ x @ 
72 0x0 +5 MEE )6:-05)) exe, pares 


2) хө 

ο ο O-BY oxo. 
| 2) 0x0 

5) охо) ©x@+2) Ox@-@E)) охо: ) @°- 
-2) 0x0 


When the sums displayed below table IV are inserted in the appropriate places, the resulting 
equations can be solved to yield 


eect NI VOVN 


ες 
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Μι _ My 
ο 521 2 τν = 38.66 
Мә м _ 12,3 Ма - 26.7 
Ly C Ly 

МЭ 197,000 

Ту 


Hence, using the value of Iy given in equations (29), 


5.21х10* Mg, = 1930 


= 
нын 
ll 


- 4.23x10° Μα = - 2670 (31) 


Ч 


Ms = ~ 1.97х10” 


To begin our discussion of these values, consider L, А compar- 
ison of the values given for Ig by equations (29) and (31) shows that 
la has been found within 3.2 percent. The nonlinear pitching moment 


Μία) = (5.21x10%)a - (4.23xK10%)a? - (1.97х107)а2 


has been superposed as the dotted curve on the true moment curve in 
figure 3. As can be seen, the agreement for this strongly nonlinear 
problem is excellent. Finally, the errors in the calculated values of 
the parameters Му and Mq are enormous, but the reason for this is well 
known and easily explained. Consider equations (13) which describe the 
motion. Eliminating q from these equations results in equation (1h); 
since the lift is linear, b, = bo = О. The constant bo, on which Mq 
and My have their principal effect, is easily interpreted physically 

as a measure of the damping in the system. The other gross aspects of— 
the response are relatively little affected by either Mg or. Μα. How- 
ever, the important quantity in bo is not Mg or Μα alone, but is their 
sum, Mg + Му. Іп other words, relatively large changes in Mg and My, 
are possible without causing any great change in the motion, Just so long 
as their sum remains constant. Thus, one may not expect to find 

Ма and Mi accurately from an experiment such as this - only their sum 
mey be relied upon. This is verified in the present exemple, since eque- 
tions (29) give the sum ман 


Mg + Mç = - ТОО 
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while equations (31) give 


Μα + Μὰ = = THO 


These two values differ by only 5.7 percent. 


It should be noted that the assumption, which has been made in 
this and the preceding example, of the linearity of the Lift is not 
necessary. If it is suspected that the lift is linear but if no def- 
inite verification of this is available, a nonlinear form such as (12) 
can be assumed for the lift, and the coefficients L, lo, La, + . ., 
can be calculated. If the lift is indeed linear, it should turn out 
that Ig, La, .. ., are small. Some limited experience has shown that 
this method does work but that the errors in the calculated parameters 
are somewhat larger than when the correct form is assumed. The reason 
for this is not known, but it appears that it may be associated with a 
tendency of the extraneous parameters I, and Lg to fit the lift curve 
to that corresponding to the original data, errors and all, at the 
expense of the smoothing operation which is necessary with this type of 
data and which is performed by the least squares process when the correct 
form is assumed. 


Example ТТТ 


We turn, finally, to a system described by a differential equation 
whose order is higher than the second. Since the higher order systems 
whose occurrence is most common appear to be those of the fourth order, 
we shall be concerned with such a system. In order to simplify the pres- 
entation, a linear system will be considered; there are no conceptual 
difficulties in the generalization to the nonlinear case. Іп addition, 
it will be assumed for simplicity that free oscillations are available 
for analysis. Thus, the system to be analyzed is assumed to be described 
by an equation of the form 


2 
atx ү ав ах + ap TH+ ау ŠZ + nox = 0 ( 32) 
at | 


dte аз at 


A solution of this equation was calculated over the interval from O 
to 2 seconds, for the following values of the coefficients 8η; 


ag = 2544.9 a, = 219.32 


132.87 аз = 2.000 


8. 2 


26 NACA TN 3288 


The result, representing the free oscillations in response to some dis- 
turbance, is tabulated in table V and presented graphically as the solid 
curve in figure 5. The sums needed for the solution of the problem are 
displayed below table V and again in table VI. The least-squares equa- 
tions for the parameters are 


24.2035 ао - 16.7225 a, - 1098.19 a» + 1603.19 a, = - 81036.4 
- 16.7225 ао + 855.538 a, + 664.495 aa - 82282.3 as = 85TTT.1 
- 1098.19 ag + 664.495 a, + 77898.7 ар - 94741.0 ag = 7533580 
1603.19 ао - 82282.3 a, - 94741.8 az + 8523670 аз = - 10567800 
solving these equations gives: - 
ёс = 3098.7 a, = 37493 
aa = 141.29 | : LE iuen gel bes ag = 3.367 


lt can be seen that these numbers are correct only to within orders of 
magnitude. On the other hand, it is not these coefficients which have 
direct. physical significance; rather, it is the damping and frequency of 
each of the components making up the oscillation which are important. 

In order to find these numbers, the following equation was set up 


A* + 3.367 AP + 111.29 AZ + 374.93 A + 3098.7 = ο 
and solved to find the roots: 


- 0,046 © 10.7 i 


- 1.64 + 4.96 i 


The true roots are obtained by solving the equation 


A* + 2.000 AS + 132.87 22 + 219.32 A+ 2544.9 = O 
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This gives 


+ 10.5 i 


- 1.00 € H.71 i 


Thus, the frequencies of the oscillation have been found quite accurately, 
as has the damping parameter of the undamped component. The only large 
effect of the errors in the coefficients ао, 81, ёд, and Әз is in the 
damping of one of the components. The apparent ill-conditioning of the 
problem with respect to this parameter is not too surprising, for after 
all, either the true or calculated value of this damping is large enough 
that the corresponding component of the motion is masked by the undamped 
component over a, good part of the run. This may be seen best, perhaps, 
from figure 5 in which the solution of equation (32), using both the 
true and calculated values of the parameters, hag been plotted. It can 
be seen that the two curves do not differ by very much, indicating that 
the fit could not be much improved. 


CONCLUDING REMARKS 


A general theory of the so-called “equations-of-motion” methods 
for the analysis of dynamical systems has been presented. It has been 
shown that, when looked at from a new point of view, all such methods 
can be generalized so as to apply to linear and nonlinear systems alike. 
Using this theory, it has also been shown how new methods can be devel- 
oped in order to satisfy the requirements of particular problems. 


One new method has been described in detail. Іп certain cases, it 
reduces to one which is very similar to the well-known. Fourier transform 
method (ref. 3) but in all cases has certain advantages over this latter 
method and over other methods heretofore used. Its superiority is based 
on two facts. First, there is the heavy dependence on the initial con- 
ditions which occurs when using most of the previously known equations- 
of-motion methods; this dependence is entirely eliminated in the new 
method. This superiority manifests itself particularly when systems of 
higher order than the second are considered. If, for definiteness, a 
fourth-order system is considered, before the Fourier transform method 
(for example) can be applied, it is necessary to evaluate the test data 
and their first three time derivatives at the initial point. Accurate 
evaluation of all these derivatives is practically impossible, however, 
with the type of data obtained from most aerodynamic experiments. 


The second fact upon which the superiority of the proposed method 
rests is that most of the equations-of-motion methods used to this time 
demand an infinitely long record for their rigorous application. For 
some years now, questions about the errors introduced into an analysis 
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of a system by the finite length of records available have been asked, - 
but answers have not been offered. The second principal advantage οἳ- 
the method described herein is that such questions are side-stepped com- 
pletely: There is no error at all from this source, since it-is assumed 
from the start that only a finite record is available. Because of this 
feature, the method avoids a further limitation of the Fourier transform 
method (apparently the most accurate of all well-known methods of this 
type), which cannot be applied at all to some systems (e.g., unstable 
ones} without time-consuming and sometimes ineffective special devices, 
since the Fourier integrals of the data simply do not exist. 


The single exception to these remarks is the derivative method 
(refs. 2 and 3). The derivative method does not-wéight the initial con- 
ditions and does not depend on an infinite interval for its application. 
In addition, the derivative method has in the past been considered as 
the only well-known method which applies to nonlinear as well as linear 
systems. (Other methods are described in references 2 and 10, but the 
derivative method appears to be the only one with such general applica- 
bility as we are discussing here.) There are, however, a number of 
very serious objections to the derivative. method. First of all, there 
is the. inordinate amount of time and labor which must be expended in its 
application, principally because of the necessity for calculating time 
derivatives of the data. Second, and most important, is the question 
of accuracy. The accurate calculation of the derivatives needed for the 
method is most-difficult, and this calculation is a large source of 
error. Besides, even if the derivatives could be computed with the 
requisite accuracy, the derivative method appears often to lead to badly 
conditioned equations, as pointed out in reference 2; because of this, 
many problems have been found for which the derivative method has been 
shown to lead to extremely large errors. The method proposed herein is 
subject to none of these weaknesses. The time required for its applica- 
tion is far less than that needed for the derivative method; in addition, 
it appears to be well suited to machine computation. Naturally, deriv- 
atives need not be calculated, and the method shares the properties of 
the Fourier transform method which cause it to lead to fairly well- 
conditioned equations. 


Ames Aeronautical Laboratory 
National Advisory Committee for Aeronautics 
Moffett Field, Calif., July 14, 1954 
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APPENDIX 
b 
NUMERICAL EVALUATION OF INTEGRALS OF THE rom Í x(t)y(t)at 
a 


In 1928, Filon (ref. 8 - see also ref. 9, pp. 67-72) published a 
generalization of Simpson's rule for evaluation of integrals of the form 


f x(t) sin wtdt 


a 


р x(t) cos wtdt 


a 


where x(t) represents numerical data. Filon'ts method, in contrast to 
Simpson's, has the distinction of giving results whose errors are inde- 
pendent of the frequency w, depending only on how closely x(t) can be 
fitted to a sequence of parabolas. This method will be generalized to 
apply to integrals of the form 


E x(t)y(t)at 


a, 


where y(t) is known exactly (for application to the method described in 
the body of this report, y(t) is one of the method functions), while 
x(t) is given tabularly. 


Suppose the interval (a,b) is divided into @h equal parts by points 
to = а <t <... < top = b, where tn} - tn = At = constant. Then, а 
formula of the form 


2h 
b 
f x(t)y(t)at а At * Tay) x(ta) (33) 
8, n=O 


will be sought, where the In are constants which depend only on the 
function y. These constants will be determined by the condition that 
the formula (33) will give the integral exactly in the cases where x(t) 
is a constant, a linear function or a quadratic function of t. Suppose 
first that the interval (a,b) is divided into two parts only by points 
to, ty, and tg; since formula (33) is to be exact if x(t) = 1, (t - ti), 


or (t - t,)%, we have 
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t 
roly) + Гу(у) + Taly) = = Í š у(%)4% 
to 


m 
- (at)roly) + (δε)ταίν) = πε Í (9 (Oe (34) 
to 
1 ба - 
(д 6) г,(у) + (А Taly) = ZZ (t - ον) y(t)at 
to 


Equations (34) can be solved for Го, Γι, Га to obtain 


tə | A To 
TERME -o£Q)2 m — - t.)y(t)at 
roly) = = L (t = ενας - аса L (5 - t,)y(4) 


tə te К | 
1 L : 2 


DA y) = 


τα - ta 
L _ +. )®у(%)а% 1 t - t,)y(t)at 
USE 4 (t - tı) y(t)ät + TTE L ( 1)y(t) 


Now, if (a,b) is divided into 2h(>2) parts, the integral is written as 
follows: | : "M 1 E E: А 


t t 
Г x(t)y(t)àt = 4 Б х(5)у(5)45 + І : x(t)y(t)dt ο ο gr 


O 2 


teh 
/ x(t)y(*)at = (36) 
Ὄργι--5 | 
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and equations (35) can be used to evaluate each of the integrals on the 


right-hand side. Calling 


t 
1 p+ 
Јо(у) = zz / y(t)at 


t 
Ko(y) - —— 1 ди (t - to)y(t)dat 


alat) TR 
БУ) = а Pte 
2 (At) 


tosa 


we obtain from (33), (35) and (36) that 


b 2h 
/᾽ κωνίοαν rat у m(y)x(to) 
8. | 1150 
where 
Го(у) = Τι. (y) ος K, (y) 
Pag) = J gp 3A! ” ӨТ-ьош1 СУ)» p = 1, 2, +. + 3 h 


Tan (y) = Lop-1 (7) + Kop1ly) T Ləp+1[ y) a Kopy1ly)»> 


p = 1, 2, . +. +> h - 1 
Tay) = Lan - (y) + Kon -11y) 


Tt should be noted that if y(t) is identically unity, 


| 
[Ὁ 


Jp(1) = 
Κρί 1) 
Lo(1) = 1/3 


i 
o 


(37) 
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where the relations to.1 = tp - At, ©р+ = tp + At have been used 
repeatedly. Hence, from equations (37), 


го(1) = 1/3 
Γαρ-ι(1) = 4/3, p = 1, 2, . . ., h 
roo (1) = 2/3, p = 1, 2, . ае 
Гор 1) = 1/3 Чин 


which exactly describes Simpson's rule. 


Equations (37), with y(t) being chosen equal to the method func- 
tions, were used to calculate the numbers displayed in tables I, III, 
and V. 
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TABLE І.- CALCULATIONS OF THE INTEGRALS NEEDED FOR EXAMPLE I - Concluded 
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TABLE V.- CALCULATION OF THE INTEGRALS NEEDED FOR EXAMPLE III - Concluded 
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Figure 1.- The nonlinear pitching-moment curve for example I. 
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Figure 2.- "Test data” for example T. 
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Figure 3.- The nonlinear pitching-moment curve for example ΤΙ, 
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Figure 4,- "Test: data” for example II, 
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Figure 5.- "Test data" for example III. 
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